data <- clean_ISSP_data

# Adjust Variable Types
data <- data %>%  mutate(high_flex_pref = as.numeric(high_flex_pref), high_meaning_pref = as.numeric(high_meaning_pref))

# Create subdata for Men & Women
data_mothers <- subset(data, woman == 1)
data_fathers <- subset(data, woman == 0)

# Function for weighted mean with CI
source(file.path(dir, "DataCode/Code/Functions/weighted_mean_ci.R"))

# Amenities by Gender
mean_flex_pref <- weighted_mean_ci(data, high_flex_pref, weight, woman)
mean_meaning_pref <- weighted_mean_ci(data, high_meaning_pref, weight, woman)

mean_pref_summary <- mean_flex_pref %>%
  rename(flex_mean = mean, flex_lb = lb, flex_ub = ub) %>%
  inner_join(mean_meaning_pref %>% rename(meaning_mean = mean, meaning_lb = lb, meaning_ub = ub), by = "woman")

# Amenities by children: Mothers
mean_flex_pref_m <- weighted_mean_ci(data_mothers, high_flex_pref, weight, children)
mean_meaning_pref_m <- weighted_mean_ci(data_mothers, high_meaning_pref, weight, children)

mean_pref_summary_m <- mean_flex_pref_m %>%
  rename(flex_mean = mean, flex_lb = lb, flex_ub = ub) %>%
  inner_join(mean_meaning_pref_m %>% rename(meaning_mean = mean, meaning_lb = lb, meaning_ub = ub), by = "children")

# Amenities by children: Fathers
mean_flex_pref_f <- weighted_mean_ci(data_fathers, high_flex_pref, weight, children)
mean_meaning_pref_f <- weighted_mean_ci(data_fathers, high_meaning_pref, weight, children)

mean_pref_summary_f <- mean_flex_pref_f %>%
  rename(flex_mean = mean, flex_lb = lb, flex_ub = ub) %>%
  inner_join(mean_meaning_pref_f %>% rename(meaning_mean = mean, meaning_lb = lb, meaning_ub = ub), by = "children")


# Prepare for Plots
plot_data_pref <- mean_pref_summary %>%
  pivot_longer(cols = c(flex_mean, meaning_mean), names_to = "variable", values_to = "mean") %>%
  mutate(
    amenity = case_when(
      variable == "meaning_mean" ~ if_else(woman == 1, 0.10, 0),
      variable == "flex_mean" ~ if_else(woman == 1, 1, 0.90)),
    lb = case_when(
      variable == "meaning_mean" ~ meaning_lb,
      variable == "flex_mean" ~ flex_lb),
    
    ub = case_when(
      variable == "meaning_mean" ~ meaning_ub,
      variable == "flex_mean" ~ flex_ub))

plot_data_pref_m <- mean_pref_summary_m %>%
  pivot_longer(cols = c(flex_mean, meaning_mean), names_to = "variable", values_to = "mean") %>%
  mutate(
    amenity = case_when(
      variable == "meaning_mean" ~ if_else(children == 1, 0.35, 0.25),
      variable == "flex_mean" ~ if_else(children == 1, 1.25, 1.15)
      ),
    lb = case_when(
      variable == "meaning_mean" ~ meaning_lb,
      variable == "flex_mean" ~ flex_lb
      ),
    ub = case_when(
      variable == "meaning_mean" ~ meaning_ub,
      variable == "flex_mean" ~ flex_ub)
    )


plot_data_pref_f <- mean_pref_summary_f %>%
  pivot_longer(cols = c(flex_mean, meaning_mean), names_to = "variable", values_to = "mean") %>%
  mutate(
    amenity = case_when(
      variable == "meaning_mean" ~ if_else(children == 1, 0.60, 0.50),
      variable == "flex_mean" ~ if_else(children == 1, 1.50, 1.40)),
    lb = case_when(
      variable == "meaning_mean" ~ meaning_lb,
      variable == "flex_mean" ~ flex_lb),
    ub = case_when(
      variable == "meaning_mean" ~ meaning_ub,
      variable == "flex_mean" ~ flex_ub))

bw_palette <- c("grey10", "grey30", "grey50", "grey70", "grey85", "grey95")

plot_pref_amenities_ISSP <- ggplot() + theme_bw() +
  geom_col(data = plot_data_pref, aes(x = amenity, y = mean, fill = c("Men", "Men", "Women", "Women")), width = 0.1, color = "black", show.legend = TRUE) +
  geom_errorbar(data = plot_data_pref, aes(x = amenity, y = mean, ymin = lb, ymax = ub), width = 0.025) +
  geom_col(data = plot_data_pref_m, aes(x = amenity, y = mean, fill = c("No Mother", "No Mother", "Mother", "Mother")), width = 0.1, color = "black", show.legend = TRUE) +
  geom_errorbar(data = plot_data_pref_m, aes(x = amenity, y = mean, ymin = lb, ymax = ub), width = 0.025) +
  geom_col(data = plot_data_pref_f, aes(x = amenity, y = mean, fill = c("No Father", "No Father", "Father", "Father")), width = 0.1, color = "black", show.legend = TRUE) +
  geom_errorbar(data = plot_data_pref_f, aes(x = amenity, y = mean, ymin = lb, ymax = ub), width = 0.025) +
  xlab("") + ylab("Percentage Highly Important") +
  scale_x_continuous(breaks = c(0.30, 1.25), labels = c("Work Meaning", "Schedule Adaptability")) +
  theme(
    legend.key = element_rect(fill = "white", colour = "white"), 
    legend.position =  "bottom", 
    legend.box = "horizontal", 
    legend.text = element_text(size = 12),
    axis.text.x = element_text(size = 14, color = "black"),  # Increase x-axis label size
    axis.text.y = element_text(size = 14, color = "black"),  # Increase y-axis label size
    axis.title.y = element_text(size = 14),  # Increase x-axis label size
    axis.ticks.x = element_blank()) + labs(color = NULL, fill = NULL) + 
  scale_fill_manual(values = bw_palette, breaks = c("Men", "Women", "No Mother", "Mother", "No Father", "Father"))  + 
  guides(color = guide_legend(override.aes = list(size = 7), name = NULL))

save_PlotAmenityPrefs <- file.path(graph_dir, "AmenityPrefHeterogeneityGenderChild_ISSP.pdf")

# Save as PNG
pdf(save_PlotAmenityPrefs) #, width = 800, height = 600)
print(plot_pref_amenities_ISSP)
dev.off()

rm("data", "data_mothers", "data_fathers", "mean_flex_pref", "mean_meaning_pref", "mean_pref_summary", "mean_flex_pref_m", "mean_meaning_pref_m",  "mean_pref_summary_m", "mean_flex_pref_f", "mean_meaning_pref_f",  "mean_pref_summary_f")
